%% This function estimate the Chen and Szroeter (2012) test
%% Inputs:
% Y=outcome, Z=instrument, D=treatment, set P=1 for step-at-unity P=2 for normal and P=3 for logistic, K=tuning parameter (possible choices sqrt(n/log(n)), sqrt(n/(2*log(log(n))))) for  des=1 if Y is descrete, cl=
% confidence level
%% Outputs
% Rc=1 if the test rejects the null
% Pc=P-value

function [Rc,Pc]=chszsim(Y,Z,D,P,K,des,cl)
%Compute the sample size
n=length(Y);
% Estimate the constraints
P11=sum(Z==1&D==1)/sum(Z==1);
P10=sum(Z==0&D==1)/sum(Z==0);
P01=1-P11;
P00=1-P10;
%Generate  the 4 conditional outcomes  Y|D=d,Z=z
y11=sort(Y(Z==1&D==1));
y01=sort(Y(Z==0&D==1));
y10=sort(Y(Z==1&D==0));
y00=sort(Y(Z==0&D==0));
%Estimate E(Y|D=1,Z=0) and E(Y|D=0,Z=1)
my01=mean(y01);
my10=mean(y10);
%Estimate the proportion of always takers in the subsample D=1, Z=1
q=(P10/P11)*(P10/P11<=1)+(P10/P11>1);
%Estimate the proportion of never takers in the subsample D=0, Z=0
r=(P01/P00)*(P01/P00<=1)+(P01/P00>1);

%Estimate the bounds for the always and never takers if bin=1 we use the correction for descrete outcomes described in Huber and Mellace (2010)
if des==0
    
    y11u=mean(y11(y11>=quantile(y11,1-q)));
    y11l=mean(y11(y11<=quantile(y11,q)));
    y00u=mean(y00(y00>=quantile(y00,1-r)));
    y00l=mean(y00(y00<=quantile(y00,r)));
else
    n11=length(y11);
    n00=length(y00);
    y11u=mean(y11(round(n11-q*n11)+(round(n11-q*n11)==0):end));
    y11l=mean(y11(1:round(q*n11)+(round(q*n11)==0)));
    
    
    y00u=mean(y00(round(n00-r*n00)+(round(n00-r*n00)==0):end));
    y00l=mean(y00(1:round(r*n00)+(round(r*n00)==0)));
end

%Estimate the moment inequality constraints
mu1=my01-y11l;
mu2=y11u-my01;
mu3=my10-y00l;
mu4=y00u-my10;


mu=[mu1;mu2;mu3;mu4];
% Estimate the asymptotic variance of sqrt(n)*(mu) (vct is a fuction see the
% related file)

V=vct(Y,Z,D,499);
% Take the main diagonal of the asymptotic variance
theta=1./sqrt(diag(V));
% 
% switch P
%     case 1 % Step-at-unity
%         phi=K*theta.*mu<=1;
%         % adjustment term
%         Lambda =-normpdf(sqrt(n)/K)*(ones(length(mu),1));
%         
%     case 2 % Normal
%         phi=1-normcdf(K*theta.*mu);
%         % adjustment term
%         Lambda =(-normpdf(K*theta.*mu)*K)/sqrt(n);
%     case 3 % Logistic
%         phi=(1+exp(K*theta.*mu)).^(-1);
%         % adjustment term
%         Lambda =((-exp(K*theta.*mu)./(1+exp(K*theta.*mu)).^2)*K)/sqrt(n);
% end
% 
% % Create a diagonal matrix with the variances of the elements of mu
% delta=diag(theta);
% % Compute the numerator of the test statistics
% Q1=sqrt(n)*phi'*delta*mu-ones(1,length(mu))*Lambda;
% % Compute the denominator of the test statistics
% Q2=sqrt(phi'*delta*V*delta*phi);
% % Compute the P-value
% if Q2==0
%     Pc=1;
% else
%     Pc=normcdf(Q1/Q2);
% end
% % Reject if Pc<cl
% Rc=Pc<cl;
mu=mu(theta<Inf & theta>-Inf & theta ~=0);
V=V(theta<Inf & theta>-Inf & theta ~=0,theta<Inf & theta>-Inf & theta ~=0);
theta=theta(theta<Inf & theta>-Inf & theta ~=0);
switch P
    case 1 % Step-at-unity
        phi=K*theta.*mu<=1;
        % adjustment term
        Lambda =-normpdf(sqrt(n)/K)*(ones(length(mu),1));
        
    case 2 % Normal
        phi=1-normcdf(K*theta.*mu);
        % adjustment term
        Lambda =(-normpdf(K*theta.*mu)*K)/sqrt(n);
    case 3 % Logistic
        phi=(1+exp(K*theta.*mu)).^(-1);
        % adjustment term
        Lambda =((-exp(K*theta.*mu)./(1+exp(K*theta.*mu)).^2))*(K/sqrt(n));
end

% Create a diagonal matrix with the variances of the elements of mu
delta=diag(theta);
% Compute the numerator of the test statistics
Q1=sqrt(n)*phi'*delta*mu-ones(1,length(mu))*Lambda;
% Compute the denominator of the test statistics
Q2=sqrt(phi'*delta*V*delta*phi);
% Compute the P-value
if Q2==0 
    Pc=1;
else
    Pc=normcdf(Q1/Q2);
end
% Reject if Pc<cl
Rc=Pc<cl;